DSi fluxes, FN concentrations, etc from WRTDS over time
1 Setup and data wrangling
reading files and adding/checking LTER names
sisyn <- read_csv("20201111_masterdata_RAW.csv")
wrtds_annual <- read_csv("WRTDS_AnnualResults_AllSites_052621.csv") #%>%
#rename(site=Site)
wrtds_trends <- read_csv("Si_EGRETCi_Trends_AllSites.csv") %>%
rename(site=Site)
sisyn_sites <- sisyn %>% select(site, LTER) %>% unique()
wrtds_annual <- left_join(wrtds_annual, sisyn_sites)
wrtds_trends <- left_join(wrtds_trends, sisyn_sites) %>%
mutate(LTER_site=paste(LTER, site, sep = "------"))
x <- wrtds_annual %>% filter(is.na(LTER))
unique(x$site)## [1] "ws1" "ws2" "ws3" "ws4" "ws5" "ws6" "ws7" "ws8" "ws9"wrtds_annual <- wrtds_annual %>%
mutate(LTER=ifelse(is.na(LTER), "HBR", LTER))
wrtds_trends <- wrtds_trends %>%
mutate(LTER=ifelse(is.na(LTER), "HBR", LTER))
z <- wrtds_annual %>% count(LTER,site)
z## # A tibble: 52 x 3
## LTER site n
## <chr> <chr> <int>
## 1 AND GSMACK 24
## 2 AND GSWS02 37
## 3 AND GSWS06 17
## 4 AND GSWS07 18
## 5 AND GSWS08 37
## 6 AND GSWS09 37
## 7 AND GSWS10 37
## 8 ARC Imnavait Weir 17
## 9 HBR ws1 48
## 10 HBR ws2 54
## # ... with 42 more rowspivoting the annual discharge/flux/conc/FNconc data and then calculating correlations and regressions across time
wrtds_long <- wrtds_annual %>%
pivot_longer(cols = c(Discharge_cms:FNFlux_106_kg_y),
names_to = "variable",
values_to = "value")
#correlations:
wrtds_corr <- wrtds_long %>%
group_by(site, LTER, variable) %>%
summarize(ri=cor(Year, value, use = "pairwise.complete.obs"), ni=n())
#regressions:
wrtds_slopes <- wrtds_long %>%
filter(!is.na(value)) %>%
group_by(site, LTER, variable) %>%
do(model = tidy(lm(value ~ Year, data = .))) %>% # only works if you include the broom::tidy in here
unnest(model)2 Correlations/slopes from annual estimates
ggplot(wrtds_corr, aes(LTER, ri)) +
stat_summary(fun = mean, geom = "point") +
stat_summary(fun.data = mean_cl_normal, geom = "pointrange", fun.args = list(mult = 1))+
geom_hline(yintercept = 0)+
xlab("LTER")+
ylab("correlation of DSi ~ year")+
stat_summary(fun.data = give.n, geom = "text", size=3, fontface="italic", position = position_nudge(x = 0.4, y=.1))+
coord_flip()+
scale_x_discrete(limits = rev)+
facet_wrap(~variable)
ggplot(wrtds_corr, aes(site, ri)) +
stat_summary(fun = mean, geom = "point") +
stat_summary(fun.data = mean_cl_normal, geom = "pointrange", fun.args = list(mult = 1))+
geom_hline(yintercept = 0)+
xlab("LTER")+
ylab("correlation of DSi ~ year")+
stat_summary(fun.data = give.n, geom = "text", size=3, fontface="italic", position = position_nudge(x = 0.4, y=.1))+
coord_flip()+
scale_x_discrete(limits = rev)+
facet_wrap(~variable)
discharge is super correlated to flux but not conc across sites and years? does this make sense?
pairs.panels(wrtds_annual)
models using correlation data
ES_wrtds <- escalc(measure="ZCOR", ri=ri, ni=ni, data=wrtds_corr)
mod1<-rma(yi=yi,vi=vi,data=filter(ES_wrtds, variable=='Discharge_cms'))
summary(mod1)##
## Random-Effects Model (k = 52; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## 1.3614 -2.7228 1.2772 5.1409 1.5272
##
## tau^2 (estimated amount of total heterogeneity): 0.0089 (SE = 0.0098)
## tau (square root of estimated tau^2 value): 0.0944
## I^2 (total heterogeneity / total variability): 17.62%
## H^2 (total variability / sampling variability): 1.21
##
## Test for Heterogeneity:
## Q(df = 51) = 63.2050, p-val = 0.1173
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.1626 0.0315 5.1603 <.0001 0.1009 0.2244 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod2<-rma(yi=yi,vi=vi,data=filter(ES_wrtds, variable=='Conc_mgL'))
summary(mod2)##
## Random-Effects Model (k = 52; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -49.4216 98.8432 102.8432 106.7069 103.0932
##
## tau^2 (estimated amount of total heterogeneity): 0.3623 (SE = 0.0809)
## tau (square root of estimated tau^2 value): 0.6019
## I^2 (total heterogeneity / total variability): 89.68%
## H^2 (total variability / sampling variability): 9.69
##
## Test for Heterogeneity:
## Q(df = 51) = 528.5065, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.0619 0.0887 0.6985 0.4849 -0.1119 0.2357
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod3<-rma(yi=yi,vi=vi,data=filter(ES_wrtds, variable=='FNConc_mgL'))
summary(mod3)##
## Random-Effects Model (k = 52; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -79.7092 159.4184 163.4184 167.2821 163.6684
##
## tau^2 (estimated amount of total heterogeneity): 1.2860 (SE = 0.2639)
## tau (square root of estimated tau^2 value): 1.1340
## I^2 (total heterogeneity / total variability): 96.86%
## H^2 (total variability / sampling variability): 31.83
##
## Test for Heterogeneity:
## Q(df = 51) = 1529.9495, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.0859 0.1601 -0.5363 0.5917 -0.3996 0.2279
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod3.1<-rma(yi=yi,vi=vi,data=filter(ES_wrtds, variable=='FNConc_mgL'),
mods = ~LTER)
summary(mod3.1)##
## Mixed-Effects Model (k = 52; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -62.3706 124.7411 146.7411 165.8555 155.5411
##
## tau^2 (estimated amount of residual heterogeneity): 1.0938 (SE = 0.2487)
## tau (square root of estimated tau^2 value): 1.0459
## I^2 (residual heterogeneity / unaccounted variability): 96.37%
## H^2 (unaccounted variability / sampling variability): 27.52
## R^2 (amount of heterogeneity accounted for): 14.94%
##
## Test for Residual Heterogeneity:
## QE(df = 42) = 1096.5069, p-val < .0001
##
## Test of Moderators (coefficients 2:10):
## QM(df = 9) = 17.5632, p-val = 0.0406
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.0112 0.4030 0.0277 0.9779 -0.7787 0.8010
## LTERARC 1.3077 1.1522 1.1349 0.2564 -0.9507 3.5660
## LTERHBR 0.3572 0.5361 0.6663 0.5052 -0.6935 1.4079
## LTERKRR(Julian) 0.0673 0.6713 0.1003 0.9201 -1.2483 1.3830
## LTERLMP(Wymore) 0.8226 1.1546 0.7124 0.4762 -1.4405 3.0856
## LTERLUQ 0.4282 0.7355 0.5822 0.5604 -1.0133 1.8698
## LTERMCM -1.3196 0.5525 -2.3885 0.0169 -2.4024 -0.2368 *
## LTERNWT 0.7491 0.7376 1.0155 0.3098 -0.6966 2.1948
## LTERSagehen(Sullivan) -0.4160 1.1385 -0.3654 0.7148 -2.6475 1.8154
## LTERUMR(Jankowski) -0.2156 0.4886 -0.4413 0.6590 -1.1733 0.7420
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 12.1 meta-regression models
wrtds_slopes <- left_join(wrtds_slopes, wrtds_corr) %>%
mutate(mi=1) %>%
filter(term=='Year')
wrtds_slopes <- wrtds_slopes %>%
arrange(LTER)
wrtds_slopes_ES_Si <- escalc(measure='ZPCOR',
ti=statistic,
ni=ni,
mi=mi,
data=filter(wrtds_slopes, variable=='Conc_mgL'))
wrtds_slopes_ES_discharge <- escalc(measure='ZPCOR',
ti=statistic,
ni=ni,
mi=mi,
data=filter(wrtds_slopes, variable=='Discharge_cms'))
wrtds_slopes_ES_flux <- escalc(measure='ZPCOR',
ti=statistic,
ni=ni,
mi=mi,
data=filter(wrtds_slopes, variable=='Flux_106_kg_y'))
mod_Si1<-rma.mv(yi,vi,
data=filter(wrtds_slopes_ES_Si),
random = ~1|site)
summary(mod_Si1)##
## Multivariate Meta-Analysis Model (k = 52; method: REML)
##
## logLik Deviance AIC BIC AICc
## -49.4102 98.8205 102.8205 106.6841 103.0705
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 0.3641 0.6034 52 no site
##
## Test for Heterogeneity:
## Q(df = 51) = 549.1464, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.0621 0.0886 0.7005 0.4836 -0.1116 0.2358
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1forest(mod_Si1)
mod_Si2<-rma.mv(yi,vi,
data=filter(wrtds_slopes_ES_Si))
summary(mod_Si2)##
## Multivariate Meta-Analysis Model (k = 52; method: REML)
##
## logLik Deviance AIC BIC AICc
## -240.9324 481.8647 483.8647 485.7965 483.9464
##
## Variance Components: none
##
## Test for Heterogeneity:
## Q(df = 51) = 549.1464, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.0523 0.0277 1.8867 0.0592 -0.0020 0.1066 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1forest(mod_Si2)
mod_discharge1<-rma.mv(yi,vi,
data=filter(wrtds_slopes_ES_discharge),
random = ~1|site)
summary(mod_discharge1)##
## Multivariate Meta-Analysis Model (k = 52; method: REML)
##
## logLik Deviance AIC BIC AICc
## 1.3103 -2.6207 1.3793 5.2430 1.6293
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 0.0109 0.1043 52 no site
##
## Test for Heterogeneity:
## Q(df = 51) = 66.1827, p-val = 0.0749
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.1628 0.0316 5.1462 <.0001 0.1008 0.2248 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1forest(mod_discharge1)
mod_flux1<-rma.mv(yi,vi,
data=filter(wrtds_slopes_ES_flux),
random = ~1|site)
summary(mod_flux1)##
## Multivariate Meta-Analysis Model (k = 52; method: REML)
##
## logLik Deviance AIC BIC AICc
## -5.8999 11.7998 15.7998 19.6635 16.0498
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 0.0335 0.1831 52 no site
##
## Test for Heterogeneity:
## Q(df = 51) = 93.2029, p-val = 0.0003
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.2005 0.0382 5.2491 <.0001 0.1256 0.2754 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1forest(mod_flux1)
viz_forest(x = wrtds_slopes_ES_Si[, c("yi", "vi")],
study_labels = wrtds_slopes_ES_Si[, "site"],
xlab = "Slope of Si by year (Z-transformed)",
variant = "thick",
group=wrtds_slopes_ES_Si[, "LTER"],
method = "REML",
annotate_CI = TRUE,
summary_label = c("Summary(LTER = AND)",
"Summary(LTER = ARC)",
"Summary(LTER = KRR(Julian))",
"Summary(LTER = HBR)",
"Summary(LTER = LMP(Wymore))",
"Summary(LTER = LUQ)",
"Summary (LTER = MCM)",
"Summary (LTER = NWT)",
"Summary(LTER = Sagehen(Sullivan))",
"Summary(LTER = UMR(Jankowski))"))
2.2 forest plots for modeled DSi conc by LTER
a <- viz_forest(x = filter(wrtds_slopes_ES_Si, LTER=="AND")[, c("yi", "vi")],
study_labels = filter(wrtds_slopes_ES_Si, LTER=="AND")[, "site"],
xlab = "ANDREWS: Slope of DSi by year (Z-transf.)",
variant = "rain")
a 
b <- viz_forest(x = filter(wrtds_slopes_ES_Si, LTER=="ARC")[, c("yi", "vi")],
study_labels = filter(wrtds_slopes_ES_Si, LTER=="ARC")[, "site"],
xlab = "ARCTIC: Slope of DSi by year (Z-transf.)",
variant = "rain")
b 
c <- viz_forest(x = filter(wrtds_slopes_ES_Si, LTER=="HBR")[, c("yi", "vi")],
study_labels = filter(wrtds_slopes_ES_Si, LTER=="HBR")[, "site"],
xlab = "HUBBARD BROOK: Slope of DSi by year (Z-transf.)",
variant = "rain")
c
d <- viz_forest(x = filter(wrtds_slopes_ES_Si, LTER=="KRR(Julian)")[, c("yi", "vi")],
study_labels = filter(wrtds_slopes_ES_Si, LTER=="KRR(Julian)")[, "site"],
xlab = "KISSIMMEE: Slope of DSi by year (Z-transf.)",
variant = "rain")
d
e <- viz_forest(x = filter(wrtds_slopes_ES_Si, LTER=="LMP(Wymore)")[, c("yi", "vi")],
study_labels = filter(wrtds_slopes_ES_Si, LTER=="LMP(Wymore)")[, "site"],
xlab = "LAMPREY: Slope of DSi by year (Z-transf.)",
variant = "rain")
e
f <- viz_forest(x = filter(wrtds_slopes_ES_Si, LTER=="LUQ")[, c("yi", "vi")],
study_labels = filter(wrtds_slopes_ES_Si, LTER=="LUQ")[, "site"],
xlab = "LUQUILLO: Slope of DSi by year (Z-transf.)",
variant = "rain")
f
g <- viz_forest(x = filter(wrtds_slopes_ES_Si, LTER=="MCM")[, c("yi", "vi")],
study_labels = filter(wrtds_slopes_ES_Si, LTER=="MCM")[, "site"],
xlab = "MCMURDO: Slope of DSi by year (Z-transf.)",
variant = "rain")
g
h <- viz_forest(x = filter(wrtds_slopes_ES_Si, LTER=="NWT")[, c("yi", "vi")],
study_labels = filter(wrtds_slopes_ES_Si, LTER=="NWT")[, "site"],
xlab = "NIWOT: Slope of DSi by year (Z-transf.)",
variant = "rain")
h
i <- viz_forest(x = filter(wrtds_slopes_ES_Si, LTER=="Sagehen(Sullivan)")[, c("yi", "vi")],
study_labels = filter(wrtds_slopes_ES_Si, LTER=="Sagehen(Sullivan)")[, "site"],
xlab = "Sagehen: Slope of DSi by year (Z-transf.)",
variant = "rain")
i
j <- viz_forest(x = filter(wrtds_slopes_ES_Si, LTER=="UMR(Jankowski)")[, c("yi", "vi")],
study_labels = filter(wrtds_slopes_ES_Si, LTER=="UMR(Jankowski)")[, "site"],
xlab = "Sagehen: Slope of DSi by year (Z-transf.)",
variant = "rain")
j
ggarrange(a,b,c,d,e,f,g,h,i,j, ncol=5, nrow=2)
viz_forest(x = wrtds_slopes_ES_Si[, c("yi", "vi")],
xlab = "Slope of DSi by year (Z-transformed)",
group=wrtds_slopes_ES_Si[, "LTER"],
variant = "rain",
study_labels = wrtds_slopes_ES_Si[, "site"],
summary_label = c("Summary(LTER = AND)",
"Summary(LTER = ARC)",
"Summary(LTER = KRR(Julian))",
"Summary(LTER = HBR)",
"Summary(LTER = LMP(Wymore))",
"Summary(LTER = LUQ)",
"Summary (LTER = MCM)",
"Summary (LTER = NWT)",
"Summary(LTER = Sagehen(Sullivan))",
"Summary(LTER = UMR(Jankowski))"))
2.3 forest plots for discharge
a <- viz_forest(x = filter(wrtds_slopes_ES_discharge, LTER=="AND")[, c("yi", "vi")],
study_labels = filter(wrtds_slopes_ES_discharge, LTER=="AND")[, "site"],
xlab = "ANDREWS: Slope of discharge by year (Z-transf.)",
variant = "rain")
a 
b <- viz_forest(x = filter(wrtds_slopes_ES_discharge, LTER=="ARC")[, c("yi", "vi")],
study_labels = filter(wrtds_slopes_ES_discharge, LTER=="ARC")[, "site"],
xlab = "ARCTIC: Slope of discharge by year (Z-transf.)",
variant = "rain")
b 
c <- viz_forest(x = filter(wrtds_slopes_ES_discharge, LTER=="HBR")[, c("yi", "vi")],
study_labels = filter(wrtds_slopes_ES_discharge, LTER=="HBR")[, "site"],
xlab = "HUBBARD BROOK: Slope of discharge by year (Z-transf.)",
variant = "rain")
c
d <- viz_forest(x = filter(wrtds_slopes_ES_discharge, LTER=="KRR(Julian)")[, c("yi", "vi")],
study_labels = filter(wrtds_slopes_ES_discharge, LTER=="KRR(Julian)")[, "site"],
xlab = "KISSIMMEE: Slope of discharge by year (Z-transf.)",
variant = "rain")
d
e <- viz_forest(x = filter(wrtds_slopes_ES_discharge, LTER=="LMP(Wymore)")[, c("yi", "vi")],
study_labels = filter(wrtds_slopes_ES_discharge, LTER=="LMP(Wymore)")[, "site"],
xlab = "LAMPREY: Slope of discharge by year (Z-transf.)",
variant = "rain")
e
f <- viz_forest(x = filter(wrtds_slopes_ES_discharge, LTER=="LUQ")[, c("yi", "vi")],
study_labels = filter(wrtds_slopes_ES_discharge, LTER=="LUQ")[, "site"],
xlab = "LUQUILLO: Slope of discharge by year (Z-transf.)",
variant = "rain")
f
g <- viz_forest(x = filter(wrtds_slopes_ES_discharge, LTER=="MCM")[, c("yi", "vi")],
study_labels = filter(wrtds_slopes_ES_discharge, LTER=="MCM")[, "site"],
xlab = "MCMURDO: Slope of discharge by year (Z-transf.)",
variant = "rain")
g
h <- viz_forest(x = filter(wrtds_slopes_ES_discharge, LTER=="NWT")[, c("yi", "vi")],
study_labels = filter(wrtds_slopes_ES_discharge, LTER=="NWT")[, "site"],
xlab = "NIWOT: Slope of discharge by year (Z-transf.)",
variant = "rain")
h
i <- viz_forest(x = filter(wrtds_slopes_ES_discharge, LTER=="Sagehen(Sullivan)")[, c("yi", "vi")],
study_labels = filter(wrtds_slopes_ES_discharge, LTER=="Sagehen(Sullivan)")[, "site"],
xlab = "Sagehen: Slope of discharge by year (Z-transf.)",
variant = "rain")
i
j <- viz_forest(x = filter(wrtds_slopes_ES_discharge, LTER=="UMR(Jankowski)")[, c("yi", "vi")],
study_labels = filter(wrtds_slopes_ES_discharge, LTER=="UMR(Jankowski)")[, "site"],
xlab = "Sagehen: Slope of discharge by year (Z-transf.)",
variant = "rain")
j
ggarrange(a,b,c,d,e,f,g,h,i,j, ncol=5, nrow=2)
viz_forest(x = wrtds_slopes_ES_discharge[, c("yi", "vi")],
xlab = "Slope of discharge by year (Z-transformed)",
group=wrtds_slopes_ES_discharge[, "LTER"],
study_labels = wrtds_slopes_ES_discharge[, "site"],
variant = "rain",
summary_label = c("Summary(LTER = AND)",
"Summary(LTER = ARC)",
"Summary(LTER = KRR(Julian))",
"Summary(LTER = HBR)",
"Summary(LTER = LMP(Wymore))",
"Summary(LTER = LUQ)",
"Summary (LTER = MCM)",
"Summary (LTER = NWT)",
"Summary(LTER = Sagehen(Sullivan))",
"Summary(LTER = UMR(Jankowski))"))
2.4 forest plots for flux
a <- viz_forest(x = filter(wrtds_slopes_ES_flux, LTER=="AND")[, c("yi", "vi")],
study_labels = filter(wrtds_slopes_ES_flux, LTER=="AND")[, "site"],
xlab = "ANDREWS: slope of DSi flux by year (Z-transf.)",
variant = "rain")
a 
b <- viz_forest(x = filter(wrtds_slopes_ES_flux, LTER=="ARC")[, c("yi", "vi")],
study_labels = filter(wrtds_slopes_ES_flux, LTER=="ARC")[, "site"],
xlab = "ARCTIC: slope of DSi flux by year (Z-transf.)",
variant = "rain")
b 
c <- viz_forest(x = filter(wrtds_slopes_ES_flux, LTER=="HBR")[, c("yi", "vi")],
study_labels = filter(wrtds_slopes_ES_flux, LTER=="HBR")[, "site"],
xlab = "HUBBARD BROOK: slope of DSi flux by year (Z-transf.)",
variant = "rain")
c
d <- viz_forest(x = filter(wrtds_slopes_ES_flux, LTER=="KRR(Julian)")[, c("yi", "vi")],
study_labels = filter(wrtds_slopes_ES_flux, LTER=="KRR(Julian)")[, "site"],
xlab = "KISSIMMEE: slope of DSi flux by year (Z-transf.)",
variant = "rain")
d
e <- viz_forest(x = filter(wrtds_slopes_ES_flux, LTER=="LMP(Wymore)")[, c("yi", "vi")],
study_labels = filter(wrtds_slopes_ES_flux, LTER=="LMP(Wymore)")[, "site"],
xlab = "LAMPREY: slope of DSi flux by year (Z-transf.)",
variant = "rain")
e
f <- viz_forest(x = filter(wrtds_slopes_ES_flux, LTER=="LUQ")[, c("yi", "vi")],
study_labels = filter(wrtds_slopes_ES_flux, LTER=="LUQ")[, "site"],
xlab = "LUQUILLO: slope of DSi flux by year (Z-transf.)",
variant = "rain")
f
g <- viz_forest(x = filter(wrtds_slopes_ES_flux, LTER=="MCM")[, c("yi", "vi")],
study_labels = filter(wrtds_slopes_ES_flux, LTER=="MCM")[, "site"],
xlab = "MCMURDO: slope of DSi flux by year (Z-transf.)",
variant = "rain")
g
h <- viz_forest(x = filter(wrtds_slopes_ES_flux, LTER=="NWT")[, c("yi", "vi")],
study_labels = filter(wrtds_slopes_ES_flux, LTER=="NWT")[, "site"],
xlab = "NIWOT: slope of DSi flux by year (Z-transf.)",
variant = "rain")
h
i <- viz_forest(x = filter(wrtds_slopes_ES_flux, LTER=="Sagehen(Sullivan)")[, c("yi", "vi")],
study_labels = filter(wrtds_slopes_ES_flux, LTER=="Sagehen(Sullivan)")[, "site"],
xlab = "Sagehen: slope of DSi flux by year (Z-transf.)",
variant = "rain")
i
j <- viz_forest(x = filter(wrtds_slopes_ES_flux, LTER=="UMR(Jankowski)")[, c("yi", "vi")],
study_labels = filter(wrtds_slopes_ES_flux, LTER=="UMR(Jankowski)")[, "site"],
xlab = "Sagehen: slope of DSi flux by year",
variant = "rain")
j
ggarrange(a,b,c,d,e,f,g,h,i,j, ncol=5, nrow=2)
viz_forest(x = wrtds_slopes_ES_flux[, c("yi", "vi")],
xlab = "slope of DSi flux by year (Z-transformed)",
group=wrtds_slopes_ES_flux[, "LTER"],
study_labels = wrtds_slopes_ES_flux[, "site"],
variant = "rain",
summary_label = c("Summary(LTER = AND)",
"Summary(LTER = ARC)",
"Summary(LTER = KRR(Julian))",
"Summary(LTER = HBR)",
"Summary(LTER = LMP(Wymore))",
"Summary(LTER = LUQ)",
"Summary (LTER = MCM)",
"Summary (LTER = NWT)",
"Summary(LTER = Sagehen(Sullivan))",
"Summary(LTER = UMR(Jankowski))"))
# fit <- nlme(data=wrtds, FNConc_mgL ~ Year, random = ~LTER)
#
# fit <- lmer(data=wrtds2, FNConc_mgL ~ Year +1|LTER)
# summary(fit)3 Using WRTDS trend estimates
hmm some have 95% CIs in the 90 million-ish range so let’s scrap those for now
ggplot(wrtds_trends, aes(LTER_site, estC)) +
geom_point()+
geom_linerange(aes(ymin = lowC95, ymax = upC95))+
coord_flip()+
geom_hline(yintercept = 0)+
scale_x_discrete(limits = rev)
wrtds_trends %>%
arrange(lowC95)## # A tibble: 52 x 33
## LTER site rejectC pValC estC lowC upC lowC50 upC50 lowC95
## <chr> <chr> <lgl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 UMR M078~ FALSE 0.146 -0.446 -9.46e+7 0.236 -8.71e-1 -0.449 -9.53e7
## 2 UMR_~ CH00~ FALSE 0.503 0.305 -4.94e-1 1.15 -2.63e-3 0.720 -1.47e6
## 3 UMR_~ CU11~ FALSE 0.271 -0.771 -2.00e+0 0.743 -1.23e+0 -0.358 -2.36e5
## 4 UMR_~ LM00~ FALSE 0.584 0.146 -7.68e-1 0.906 -1.24e-1 0.536 -1.67e5
## 5 LUQ RI FALSE 0.886 1.33 -2.28e+4 5.66 -1.69e+4 2.16 -2.55e4
## 6 Sage~ Sage~ TRUE 0.0198 -1.41 -2.14e+4 -0.720 -2.09e+4 -1.21 -2.18e4
## 7 NWT MART~ TRUE 0.0944 0.263 3.73e-2 0.414 1.65e-1 0.325 -1.77e4
## 8 ARC Imna~ FALSE 0.795 0.553 -6.63e+0 13.1 -5.95e-1 0.718 -3.69e1
## 9 NWT SADD~ FALSE 0.374 5.34 -5.51e+0 10.2 1.23e+0 6.97 -1.75e1
## 10 LUQ MPR FALSE 0.834 -0.246 -4.98e+0 3.25 -2.94e+0 0.786 -6.74e0
## # ... with 42 more rows, and 23 more variables: upC95 <dbl>, likeCUp <dbl>,
## # likeCDown <dbl>, rejectF <lgl>, pValF <dbl>, estF <dbl>, lowF <dbl>,
## # upF <dbl>, lowF50 <dbl>, upF50 <dbl>, lowF95 <dbl>, upF95 <dbl>,
## # likeFUp <dbl>, likeFDown <dbl>, baseConc <dbl>, baseFlux <dbl>,
## # iBoot <dbl>, startSeed <dbl>, nBootGood <dbl>, Year.1 <dbl>, Year.2 <dbl>,
## # duration <dbl>, LTER_site <chr>wrtds_trends %>%
arrange(desc(upC95))## # A tibble: 52 x 33
## LTER site rejectC pValC estC lowC upC lowC50 upC50 lowC95
## <chr> <chr> <lgl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ARC Imna~ FALSE 0.795 0.553 -6.63e+0 13.1 -5.95e-1 0.718 -3.69e+1
## 2 NWT SADD~ FALSE 0.374 5.34 -5.51e+0 10.2 1.23e+0 6.97 -1.75e+1
## 3 LUQ RI FALSE 0.886 1.33 -2.28e+4 5.66 -1.69e+4 2.16 -2.55e+4
## 4 UMR_~ WP02~ FALSE 0.433 -1.52 -2.10e+0 1.44 -1.63e+0 -0.157 -2.37e+0
## 5 LUQ MPR FALSE 0.834 -0.246 -4.98e+0 3.25 -2.94e+0 0.786 -6.74e+0
## 6 MCM Ande~ FALSE 0.339 -0.309 -6.11e-1 2.30 -4.52e-1 -0.120 -7.45e-1
## 7 UMR_~ BK01~ FALSE 0.138 0.482 -9.85e-2 1.69 2.74e-1 1.04 -3.44e-1
## 8 KRR S65E FALSE 0.602 0.280 -2.36e+0 1.52 -2.62e-1 1.01 -3.18e+0
## 9 UMR_~ CN00~ FALSE 0.504 0.223 -4.01e-1 1.57 -1.16e-2 0.970 -6.30e-1
## 10 LUQ QS TRUE 0.0392 0.990 4.70e-1 1.64 7.72e-1 1.26 2.13e-1
## # ... with 42 more rows, and 23 more variables: upC95 <dbl>, likeCUp <dbl>,
## # likeCDown <dbl>, rejectF <lgl>, pValF <dbl>, estF <dbl>, lowF <dbl>,
## # upF <dbl>, lowF50 <dbl>, upF50 <dbl>, lowF95 <dbl>, upF95 <dbl>,
## # likeFUp <dbl>, likeFDown <dbl>, baseConc <dbl>, baseFlux <dbl>,
## # iBoot <dbl>, startSeed <dbl>, nBootGood <dbl>, Year.1 <dbl>, Year.2 <dbl>,
## # duration <dbl>, LTER_site <chr>ggplot(wrtds_trends, aes(estC))+
geom_histogram()
ggplot(wrtds_trends, aes(lowC95))+
geom_histogram()
ggplot(wrtds_trends, aes(upC95))+
geom_histogram()
wrtds_trends %>%
filter(lowC95 < -2e7)## # A tibble: 1 x 33
## LTER site rejectC pValC estC lowC upC lowC50 upC50 lowC95 upC95
## <chr> <chr> <lgl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 UMR M078~ FALSE 0.146 -0.446 -9.46e7 0.236 -0.871 -0.449 -9.53e7 0.906
## # ... with 22 more variables: likeCUp <dbl>, likeCDown <dbl>, rejectF <lgl>,
## # pValF <dbl>, estF <dbl>, lowF <dbl>, upF <dbl>, lowF50 <dbl>, upF50 <dbl>,
## # lowF95 <dbl>, upF95 <dbl>, likeFUp <dbl>, likeFDown <dbl>, baseConc <dbl>,
## # baseFlux <dbl>, iBoot <dbl>, startSeed <dbl>, nBootGood <dbl>,
## # Year.1 <dbl>, Year.2 <dbl>, duration <dbl>, LTER_site <chr>wrtds_trends %>%
filter(upC95 > 2e7)## # A tibble: 1 x 33
## LTER site rejectC pValC estC lowC upC lowC50 upC50 lowC95 upC95 likeCUp
## <chr> <chr> <lgl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ARC Imna~ FALSE 0.795 0.553 -6.63 13.1 -0.595 0.718 -36.9 5.17e7 NA
## # ... with 21 more variables: likeCDown <dbl>, rejectF <lgl>, pValF <dbl>,
## # estF <dbl>, lowF <dbl>, upF <dbl>, lowF50 <dbl>, upF50 <dbl>, lowF95 <dbl>,
## # upF95 <dbl>, likeFUp <dbl>, likeFDown <dbl>, baseConc <dbl>,
## # baseFlux <dbl>, iBoot <dbl>, startSeed <dbl>, nBootGood <dbl>,
## # Year.1 <dbl>, Year.2 <dbl>, duration <dbl>, LTER_site <chr>wrtds_trends_trimmed <- wrtds_trends %>%
filter(lowC95>-50) %>%
filter(upC95<50)summary(wrtds_trends_trimmed)## LTER site rejectC pValC
## Length:44 Length:44 Mode :logical Min. :0.01980
## Class :character Class :character FALSE:33 1st Qu.:0.09864
## Mode :character Mode :character TRUE :11 Median :0.23482
## Mean :0.33119
## 3rd Qu.:0.51405
## Max. :0.97723
##
## estC lowC upC lowC50
## Min. :-1.5187 Min. :-5.5108 Min. :-0.01031 Min. :-2.94079
## 1st Qu.:-0.1363 1st Qu.:-0.5729 1st Qu.: 0.35689 1st Qu.:-0.25621
## Median : 0.1965 Median :-0.2876 Median : 0.86837 Median :-0.05223
## Mean : 0.2363 Mean :-0.6132 Mean : 1.10630 Mean :-0.08456
## 3rd Qu.: 0.4931 3rd Qu.:-0.1114 3rd Qu.: 1.38779 3rd Qu.: 0.30061
## Max. : 5.3402 Max. : 0.9813 Max. :10.16405 Max. : 1.23391
##
## upC50 lowC95 upC95 likeCUp
## Min. :-0.58212 Min. :-17.5133 Min. : 0.02856 Min. :0.02778
## 1st Qu.:-0.01424 1st Qu.: -0.9011 1st Qu.: 0.43535 1st Qu.:0.25490
## Median : 0.42081 Median : -0.4742 Median : 1.13162 Median :0.71569
## Mean : 0.55089 Mean : -1.1159 Mean : 1.37864 Mean :0.60333
## 3rd Qu.: 0.85097 3rd Qu.: -0.2064 3rd Qu.: 1.63045 3rd Qu.:0.94059
## Max. : 6.97374 Max. : 0.8973 Max. :10.41049 Max. :0.99505
## NA's :1
## likeCDown rejectF pValF estF
## Min. :0.004951 Mode :logical Min. :0.0198 Min. :-10352.558
## 1st Qu.:0.059406 FALSE:34 1st Qu.:0.1570 1st Qu.: -0.056
## Median :0.284314 TRUE :10 Median :0.4122 Median : 0.033
## Mean :0.396673 Mean :0.4393 Mean : -303.418
## 3rd Qu.:0.745098 3rd Qu.:0.6886 3rd Qu.: 6.917
## Max. :0.972222 Max. :0.9244 Max. : 278.731
## NA's :1
## lowF upF lowF50
## Min. :-20645.685 Min. : 0.002 Min. :-13299.324
## 1st Qu.: -16.438 1st Qu.: 0.133 1st Qu.: -0.797
## Median : -0.205 Median : 0.769 Median : -0.002
## Mean : -697.895 Mean : 302.531 Mean : -432.513
## 3rd Qu.: -0.002 3rd Qu.: 53.948 3rd Qu.: 0.241
## Max. : 33.868 Max. :5488.965 Max. : 196.891
##
## upF50 lowF95 upF95 likeFUp
## Min. :-4964.677 Min. :-24081.227 Min. : 0.002 Min. :0.09406
## 1st Qu.: 0.014 1st Qu.: -37.228 1st Qu.: 0.176 1st Qu.:0.34259
## Median : 0.305 Median : -0.325 Median : 0.915 Median :0.69444
## Mean : -111.685 Mean : -816.794 Mean : 738.278 Mean :0.62305
## 3rd Qu.: 17.943 3rd Qu.: -0.012 3rd Qu.: 73.833 3rd Qu.:0.92574
## Max. : 402.210 Max. : 4.452 Max. :18713.037 Max. :0.99505
## NA's :3
## likeFDown baseConc baseFlux iBoot
## Min. :0.004951 Min. : 0.661 Min. : 0.03 Min. : 50.00
## 1st Qu.:0.074257 1st Qu.: 2.534 1st Qu.: 0.50 1st Qu.: 50.00
## Median :0.305556 Median : 4.725 Median : 3.16 Median : 50.00
## Mean :0.376951 Mean : 5.004 Mean : 1717.12 Mean : 61.93
## 3rd Qu.:0.657407 3rd Qu.: 5.935 3rd Qu.: 204.11 3rd Qu.: 59.50
## Max. :0.905941 Max. :20.101 Max. :42361.98 Max. :100.00
## NA's :3
## startSeed nBootGood Year.1 Year.2
## Min. :494817 Min. : 50.00 Min. :1964 Min. :2012
## 1st Qu.:494817 1st Qu.: 50.00 1st Qu.:1984 1st Qu.:2017
## Median :494817 Median : 50.00 Median :1995 Median :2018
## Mean :494817 Mean : 61.93 Mean :1991 Mean :2018
## 3rd Qu.:494817 3rd Qu.: 59.50 3rd Qu.:1997 3rd Qu.:2019
## Max. :494817 Max. :100.00 Max. :2003 Max. :2019
##
## duration LTER_site
## Min. :15.00 Length:44
## 1st Qu.:21.75 Class :character
## Median :22.00 Mode :character
## Mean :27.00
## 3rd Qu.:31.50
## Max. :53.00
## ggplot(wrtds_trends_trimmed, aes(site, estC)) +
geom_point()+
geom_linerange(aes(ymin = lowC95, ymax = upC95))+
coord_flip()+
geom_hline(yintercept = 0)+
scale_x_discrete(limits = rev)+
facet_wrap(~LTER, scales='free')
#trying to make summary points per LTER didn't work this way
# ggplot(wrtds_trends_trimmed) +
# geom_point(aes(site, estC))+
# geom_linerange(aes(site, ymin = lowC95, ymax = upC95))+
# coord_flip()+
# geom_hline(yintercept = 0)+
# scale_x_discrete(limits = rev)+
# stat_summary(aes(LTER, fun = mean, geom = "point")) +
# stat_summary(aes(LTER, fun.data = mean_cl_normal,
# geom = "pointrange", fun.args = list(mult = 1)))+
# facet_wrap(~LTER, scales='free')3.1 new attempt at weighted forest plots
# data stright from wrtds output
wrtds_trends_trimmed%>%
ggplot(aes(fill = LTER, x = paste(LTER, site, sep = "::"), y = estC, ymin = lowC95, ymax=upC95))+
geom_pointrange(color="black", shape = 21, size = 2)+
scale_x_discrete(limits = rev)+
coord_flip()+
geom_hline(yintercept = 0, linetype="dashed")+
theme(axis.text=element_text(size=16))+
# scale_color_manual(values = jever2color)+
# scale_fill_manual(values = jever2color)+
#scale_size_manual(values=c(1,1.5))+
#scale_shape_manual(values = c(21, 22))+
xlab("site")+
ylab("Slope of annual seasonality change")
#make SE data
wrtds_trends_trimmed <- wrtds_trends_trimmed %>%
mutate(SE_est = ((upC95-lowC95)/2)/sqrt(duration))
# do meta-analysis by LTER
mafits_trends <- wrtds_trends_trimmed %>%
group_by(LTER) %>% #group by type of nutrient manipulation and broad response
do(model = tidy(rma(estC, SE_est, data = .), conf.int=TRUE)) %>% #do weighted meta-analysis for each grouping
unnest(model)
#plot number observations and weighted data from above
num_obs <- wrtds_trends_trimmed %>%
group_by(LTER) %>%
summarize(num_obs=length(estC), mean_estimate=mean(estC))
mafits_trends## # A tibble: 9 x 9
## LTER term type estimate std.error statistic p.value conf.low conf.high
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 AND overa~ summa~ 0.0660 0.104 0.634 0.526 -0.138 0.270
## 2 HBR overa~ summa~ 0.484 0.187 2.59 0.00973 0.117 0.850
## 3 KRR overa~ summa~ 0.165 0.249 0.660 0.509 -0.324 0.654
## 4 LMP overa~ summa~ 0.851 0.512 1.66 0.0964 -0.152 1.85
## 5 LUQ overa~ summa~ 0.708 0.518 1.37 0.172 -0.308 1.72
## 6 MCM overa~ summa~ -0.113 0.0971 -1.16 0.244 -0.303 0.0773
## 7 NWT overa~ summa~ 2.50 2.52 0.991 0.322 -2.44 7.44
## 8 UMR overa~ summa~ -0.0829 0.357 -0.232 0.816 -0.783 0.617
## 9 UMR_trib overa~ summa~ 0.166 0.265 0.626 0.531 -0.354 0.686num_obs## # A tibble: 9 x 3
## LTER num_obs mean_estimate
## <chr> <int> <dbl>
## 1 AND 7 0.0419
## 2 HBR 9 0.477
## 3 KRR 4 0.167
## 4 LMP 1 0.851
## 5 LUQ 2 0.372
## 6 MCM 8 -0.122
## 7 NWT 2 2.80
## 8 UMR 5 -0.164
## 9 UMR_trib 6 -0.0415mafits_trends <- left_join(mafits_trends, num_obs)
mafits_trends%>%
ggplot(aes(x = LTER, y = estimate, ymin = conf.low, ymax=conf.high))+
geom_pointrange(color="black", shape = 21, size = 1.8, fill = "#218D3A")+
scale_x_discrete(limits = rev)+
coord_flip()+
geom_hline(yintercept = 0, linetype="dashed")+
geom_text(size=4, fontface="italic", position = position_nudge(x = 0.35, y=.25),
aes(label = num_obs))+
theme(axis.text=element_text(size=16))+
scale_color_manual(values = jever2color)+
scale_fill_manual(values = jever2color)+
#facet_wrap(~variable)+
#scale_size_manual(values=c(1,1.5))+
#scale_shape_manual(values = c(21, 22))+
xlab("LTER")+
ylab("Slope DSi change over time")
#ggsave("wrtds_trends_weighted.png")3.2 old attempts, maybe not useful
calculating a rough estimate of standard error from the 95% CIs for now, but need to make sure this is ok to do (should still give a relatively good weighting for meta-regression stuff for now/makes sense when compared to the CIs)
wrtds_trends_trimmed <- wrtds_trends_trimmed %>%
mutate(SE_est = ((upC95-lowC95)/2)/sqrt(duration))
ggplot(wrtds_trends_trimmed, aes(site, SE_est))+
geom_point()+
coord_flip()
trying to do a meta-regression model for each LTER because that’s the only way I know how to get overall weighted effect sizes for each LTER, but metafor doesn’t like whatever I’m doing
likelihood vs pval stuff to check out how those are all related
wrtds_trends %>% count(rejectC)## # A tibble: 2 x 2
## rejectC n
## <lgl> <int>
## 1 FALSE 39
## 2 TRUE 13ggplot(wrtds_trends, aes(pValC, rejectC))+
geom_point()+
geom_vline(xintercept=.05)
ggplot(wrtds_trends, aes(likeCDown, rejectC))+
geom_point()+
geom_vline(xintercept=.05)
ggplot(wrtds_trends, aes(likeCUp, rejectC))+
geom_point()+
geom_vline(xintercept=.95)
4 maps and latitude (non)effects
importing lat-long data and having a quick look- looks like Sagehen needs to swap hemispheres
latlongs <- read_csv("LongTermWatersheds_LatLong.csv") %>%
rename(site=Stream.Site)
leaflet(latlongs) %>%
addTiles() %>%
addMarkers(lng = ~Longitude, lat = ~Latitude, popup = ~ LTER)latlongs <- latlongs %>%
mutate(Longitude = ifelse(LTER=='Sagehen', -Longitude, Longitude))
leaflet(latlongs) %>%
addTiles() %>%
addMarkers(lng = ~Longitude, lat = ~Latitude, popup = ~ LTER)a <- unique(wrtds_trends$site) #usi this to check where there are differences in the colnames
b <- unique(latlongs$site)
setdiff(b,a)## [1] "TW Weir" "Toolik Inlet" "Q1" "Q2" "Q3"
## [6] "QP"b %in% a## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE TRUE TRUE
## [13] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [25] TRUE TRUE TRUE FALSE FALSE FALSE FALSE TRUE TRUE TRUE TRUE TRUE
## [37] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [49] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUEwhich(a %in% b)## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
## [51] 51 52wrtds_trends2 <- left_join(wrtds_trends, latlongs, by='site')
wrtds_slopes_ES_Si2 <- left_join(wrtds_slopes_ES_Si, latlongs, by="site")
wrtds_slopes_ES_discharge2 <- left_join(wrtds_slopes_ES_discharge, latlongs, by="site")
wrtds_slopes_ES_flux2 <- left_join(wrtds_slopes_ES_flux, latlongs, by="site")joining with slope data- here these represent untransformed values of the slopes (“real” slope units of change per time from WRTDS trend estimates I think) where all above forest plots were Z-transformed
pal <- colorBin(palette = "RdBu", bins = c(-6,-3,-2,-1,0,1,2,3,6), domain = wrtds_trends$estC, reverse = TRUE)
#pal <- colorNumeric("RdBu", wrtds_trends$estC)
leaflet(wrtds_trends2) %>%
addTiles() %>%
addCircleMarkers(lng = ~Longitude,
lat = ~Latitude,
#popup = ~ LTER.x,
color="black",
stroke=TRUE,
fillColor = ~pal(estC),
fillOpacity = 10,
opacity = .5,
label = ~paste0(wrtds_trends2$LTER.x,"::", wrtds_trends2$site, ", slope: ", estC),
radius = 10) # I don't think the legend is super helpful..
# addLegend(position = "bottomleft",
# pal = pal,
# values = range(wrtds_trends$estC),
# title = "WRTDS slope estimate (DSi~year)")4.1 meta-regression including absolute latitude as variable
plotting absolute lat vs slopes for DSi/discharge/flux
ggplot(wrtds_trends2, aes(abs(Latitude), estC))+
geom_point()+
geom_smooth(method='lm')
ggplot(wrtds_slopes_ES_Si2, aes(abs(Latitude), estimate))+
geom_point()+
geom_smooth(method='lm')+
ylab('slope DSi from annual estimates')
ggplot(wrtds_slopes_ES_discharge2, aes(abs(Latitude), estimate))+
geom_point()+
geom_smooth(method='lm')+
ylab('slope discharge from annual estimates')
ggplot(wrtds_slopes_ES_flux2, aes(abs(Latitude), estimate))+
geom_point()+
geom_smooth(method='lm')+
ylab('slope DSi flux from annual estimates')
meta-regression stuff also shows no effect of lat, but here’s where we could also plug in additional variables to see if anything modifies slopes
mod_Si3<-rma.mv(yi,vi,
data=wrtds_slopes_ES_Si2,
random = ~1|site,
mods = ~abs(Latitude))
summary(mod_Si3)##
## Multivariate Meta-Analysis Model (k = 52; method: REML)
##
## logLik Deviance AIC BIC AICc
## -47.0893 94.1786 100.1786 105.9147 100.7004
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 0.3424 0.5851 52 no site
##
## Test for Residual Heterogeneity:
## QE(df = 50) = 512.2841, p-val < .0001
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 3.8262, p-val = 0.0505
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.5588 0.2682 2.0834 0.0372 0.0331 1.0846 *
## abs(Latitude) -0.0108 0.0055 -1.9561 0.0505 -0.0217 0.0000 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod_discharge2<-rma.mv(yi,vi,
data=wrtds_slopes_ES_discharge2,
random = ~1|site,
mods = ~abs(Latitude))
summary(mod_discharge2)##
## Multivariate Meta-Analysis Model (k = 52; method: REML)
##
## logLik Deviance AIC BIC AICc
## 1.5312 -3.0624 2.9376 8.6737 3.4593
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 0.0105 0.1023 52 no site
##
## Test for Residual Heterogeneity:
## QE(df = 50) = 64.3719, p-val = 0.0832
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 1.4759, p-val = 0.2244
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.2798 0.1014 2.7599 0.0058 0.0811 0.4785 **
## abs(Latitude) -0.0026 0.0021 -1.2149 0.2244 -0.0067 0.0016
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod_flux2<-rma.mv(yi,vi,
data=wrtds_slopes_ES_flux2,
random = ~1|site,
mods = ~abs(Latitude))
summary(mod_flux2)##
## Multivariate Meta-Analysis Model (k = 52; method: REML)
##
## logLik Deviance AIC BIC AICc
## -4.9865 9.9729 15.9729 21.7090 16.4947
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2 0.0316 0.1777 52 no site
##
## Test for Residual Heterogeneity:
## QE(df = 50) = 88.8180, p-val = 0.0006
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 2.5593, p-val = 0.1096
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.3818 0.1195 3.1959 0.0014 0.1476 0.6159 **
## abs(Latitude) -0.0040 0.0025 -1.5998 0.1096 -0.0088 0.0009
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1